Molecular Line Profiles from a Core Forming in a Turbulent Cloud 



Jeong-Eun Lee 

Department of Astronomy and Space Science, Astrophysical Research Center for the 
Structure and Evolution of the Cosmos, Sejong University, 98 Cunja-Dong, Cwangjin-Cu, 

Seoul 143-747, Korea 

j eleeSse j ong . ac . kr 

and 
Jongsoo Kim 

Korea Astronomy and Space Science Institute, Hwaam-Dong, Yuseong-Gu, Daejeon 

305-348, South Korea 

j skimSkasi . re . kr 
ABSTRACT 

We calculate the evolution of molecular line profiles of HCO^ and C^^O to- 
ward a dense core that is forming inside a magnetized turbulent molecular cloud. 
Features of the profiles can be affected more significantly by coupled velocity and 
abundance structures in the outer region than those in the inner dense part of 
the core. The velocity structure at large radii is dominated by a turbulent flow 
nearby and accretion shocks onto the core, which resulting in the variation be- 
tween inward and outward motions during the evolution of the core. The chemical 
abundance structure is significantly affected by the depletion of molecules in the 
central region with high density and low temperature. During the evolution of 
the core, the asymmetry of line profiles easily changes from blue to red, and vice 
versa. According to our study, the observed reversed (red) asymmetry toward 
some starless cores could be interpreted as an intrinsic result of outward motion 
in the outer region of a dense core, which is embedded in a turbulent environment 
and still grows in density at the center. 



Subject headings: astrochemistry — ISM: molecules — ISM: clouds — stars: 
formation — MHD — methods: numerical 



- 2 - 



1. Introduction 

Stars form inside very dense molecular clouds by gravitational collapse. Starless dense 
cores, which have been considered as potential sites of future star formation, have been 
used to study the initial conditions of star formation. There is, however, no straightforward 
method to determine whether specific cores will actually form stars. Observations tell us the 
mass or density profile of a core, which provides clues to whether the core is gravitationally 
bound or unstable to collapse to form a star in the future. 

Walker et al. (1986) and Zhou et al. (1993) interpreted the asymmetry of hne profiles 
with blue stronger peaks {blue asymmetry), which have been seen in optically thick self- 
absorbed lines, as an infall signature. Gregersen & Evans (2000) and Lee et al. (2001) 
have surveyed dense starless cores with molecular fines of HCO+ 3—2 and CS 2—1. The 
infall signature is seen in about 35% and 50% of the observed starless cores of each survey, 
respectively. The infall signature, however, occurs in much larger scales across the face of 
a core opposed to the collapsing center expected from the inside-out collapse model (Shu 
1977), indicative of an overall inward motion. 

This overall inward motion in starless dense cores could be considered as a result of a 
supercritical (or slightly subcritical) ambipolar diffusion model in Ciolek & Basu (2000). The 
inward velocity profile in L1544, the canonical starless dense core, is able to be reproduced 
reasonably well with the model (Caselli et al. 2002). However, in addition to the blue 
asymmetry, the reversed asymmetry, which has a red stronger peak in the self-absorption 
feature, is also seen toward many starless cores. This reversed asymmetry cannot be caused 
by outfiows, which can be developed only after an accretion disk forms around a young stellar 
object. However, by definition, a starless core can not harbor the young stellar object. Other 
studies claim that the reversed asymmetry might be caused by rotation (Redman et al. 2004; 
Lee et al. 2007) or osciUation (Keto et al. 2006; Maret et al. 2007; Aguti et al. 2007) of 
dense cores. 

As described above, the interpretation of the velocity profiles toward starless dense cores 
is still controversial since the formation process of dense cores is not well understood. Here, 
we suggest a possible explanation of the reversed red asymmetry observed toward some 
starless cores with the theory based on turbulence since it involves not only inward but also 
outward motions in a core formed from the density fluctuations of a larger turbulent cloud 
(see, e.g., Mac Low & Kessen 2004). In addition to the dynamical process, the chemical 
process in a dense core can also affect molecular fine profiles (Rawfings & Yates 2001, Lee et 
al. 2003, 2004, 2005). Therefore, observed line profiles must be interpreted based on both 
the dynamics and chemistry of a core. 
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In this Letter, we combine dynamics and chemistry in order to study the evolution of 
molecular line profiles toward a dense core forming within a large turbulent cloud. We show 
that the velocity structure at the outer part of the core can significantly affect molecular line 
profiles. Indeed, simulated line profiles exhibit a mixture of infall and outfiow signatures as 
a function of time. Thus a snap shot of core conditions of a single object cannot provide 
a full accounting of core formation and evolution toward collapse. This Letter is organized 
as follow. The dynamical and chemical models of a core are presented in §2, results, and 
summary and discussion are given in §3 and §4, respectively. 



2. Dynamical and Chemical models of a Core 

2.1. Dynamical Evolution of a Core 

What we are interested in this study is to follow up the evolution of molecular line 
profiles of a core forming in a turbulent molecular cloud. For this purpose, we perform a 
couple of three-dimensional numerical experiments of a self-gravitating, isothermal, turbu- 
lent, magnetized molecular cloud. In fact, we integrate, as a function of time, the isothermal 
MHD and Poisson equations. The number of cells used in the experiments is 512^. A typical 
core is covered with more than 32^ cells. 

There are two parameters in the numerical experiments: one is the plasma (3, which is 
the ratio of gas pressure to magnetic pressure, and the other is the Jeans number J, which 
is the ratio of the one dimensional size of a computational box L, to the Jeans length L j. A 
mass-to-fiux ratio normalized with its critical value ^, can be written in terms of these two 
parameters = t^J^ 13/2, where we use the critical mass-to-fiux ratio (M/0)cr = (47r^G)~^/^ 
in Nagano & Nagamura (1978). Initially a uniform molecular cloud is threaded with a 
uniform magnetic field. In order to generate a turbulent fiow and keep the level of the 
turbulence inside the molecular cloud, we add velocity perturbations generated in the Fourier 
space with wavelengths that span from half to full size of the box. We adjust the input rate 
of the kinetic energy so that the root-mean-square sonic Mach number Mg of a turbulent 
fiow becomes 10 (Mac Low 1999; Stone et al. 1998). Once a fully saturated turbulent fiow 
is generated, we turn on the self-gravity of gas and follow up the evolution of a densest core 
formed in the cloud. If /3 = 0.1 and J = 4, then the n value is 2.8. So the initial molecular 
cloud is in a highly supersonic (M, >> 1), Jeans unstable (J > 1), and magnetically super- 
critical {jji> 1) state. 

Since the combined isothermal MHD and Poisson equations can be written in dimen- 
sionless form with the two parameters, we can select any units of, for example, length, time. 
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and mass in such a way that they are scalable with each other. However there is one more ad- 
ditional parameter M^, which is used to fix these units. If Larson's relation between the sizes 
and velocity dispersions of clouds (Larson 1981) holds for our model cloud, then the Mach 
number, which is the velocity dispersion of our model cloud, can determine its corresponding 
size. In fact, a Larson-type relation of the form Mg = 5(L/lpc)^/^, if we take 0.5 as the 
power index (for example, Myers 1983), is used. Then density and magnetic field strength 
could be determined from relations, n — 500( J/(L/lpc))^ cm"^, and B — 0.205(n//9)^/^ /xG 
(see also equation (7) and equation (8) in Vazquez- Semadeni et al. 2005). So if we choose 
/3 = 0.1, J — 4:, and M, = 10, the size of our computational box, the initial number density, 
and the initial magnetic field strength become 4 pc, 500 cm~^, and 14.5 //G, respectively. 

Our dynamical model has a few limitations such as turning on self-gravity after gen- 
erating a fully saturated turbulent flow, driving turbulent flows by adding velocity fields 
generated in the Fourier space, not taking into account ambipolar diffusion, and assuming 
the isothermal condition. However, they are common in numerical models of a turbulent 
molecular cloud (see, for example, Klessen et al. 2000; Ostriker et al. 2001; Li et al. 2004; 
Vazquez-Semadeni et al. 2005). Due to the limited space, we only mention the effects of the 
latter two limitations on our results. If ambipolar diffusion is included in our model, then 
the ambipolar diffusion process is more important in the central part of a core where the 
degree of ionization is lower than that of the outer part of the core. Since our synthesized 
line profiles are more sensitive to the conditions at the outer layer (see section 3), our main 
results are not likely to be significantly affected by ambipolar diffusion. We also take the 
isothermal approximation which enables us not to solve the energy equation with cooling and 
heating processes. To take into account the cooling and heating processes of gas realistically 
entails a heavy calculation because molecular cooling lines are usually optically thick. When 
line profiles arc, however, synthesized, we calculate two cases with a constant temperature 
and temperature variations of gas and dust inside the core (sec section 3). We present the 
latter case only because the line features from both cases are not much different from each 
other. 



2.2. Chemical Evolution and the HCO"*" Line Evolution 

In order to calculate the chemical evolution in the core picked-up from the dynamical 
simulation as described in the previous subsection, we adopt the evolutionary chemical model 
developed by Lee et al. (2004), in particular, the model before collapse begins. In the 
prestellar phase where this paper focuses on, the model calculates the chemical evolution at 
each grid point as density grows with time. The chemical network in the model has been 
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updated to include more recent results on the binding energy of N2 (Oberg et al. 2005) and 
the photodesorption yield of CO (Oberg et al. 2007). For the initial molecular abundances of 
the coupled evolution between dynamics and chemistry, we calculate the chemical evolution 
in a constant density of 5 x 10'^ cm^^ with the same initial atomic abundances as used in 
Lee et al. (2004) for t = 5 x 10''"^ years. For the chemical calculation, the dense core is 
assumed to be surrounded with material of Ay = 1.0 mag, where CO is self-shielded. For 
consistency, adopting the methods used in Lee et al. (2004), we also calculate the dust and 
gas temperatures with the interstellar radiation field attenuated by Ay — 1.0 mag. The time 
step for data dump of the dynamical simulation is 4 x 10^ years, during which the chemistry 
is calculated based on constant density and temperatures given at a specific dump time, and 
the result of the chemical evolution is used as initial conditions for the next dump time. 

We use the Monte-Carlo method (Choi et al. 1995) to perform the radiative transfer 
calculation for the HCO+ and C^^O hues. The collision rates for HCO"*" and CO are adopted 
from Flower (1999) and Flower & Launay (1985). For the simulation of the line profiles, 
wc assume the distance to the model core of 250 pc and the beam sizes for the Caltech 
Submillimeter Observatory 10.4 m telescope. We also include a constant microturbulent 
velocity dispersion of 0.2 km s~^. We model the HCO"'" 3—2 and 4—3 transitions, which 
have been used to study infall motion in low mass cores with and without young stellar 
objects (Gregersen et al. 1997; Gregersen & Evans 2000). The optically thin C^^O 3-2 hne 
is also modeled to be compared to the HCO+ line profiles. 



3. Results 

We use a one-dimensional (radial) chemical network coupled with a radiative transfer 
model to simulate line profiles. For this purpose we have to reduce the three-dimensional 
density and velocity fields into one-dimensional radial distributions using the following pro- 
cedures. First, we find the position of a density peak of the whole computational domain, 
which is taken as the origin of the radial coordinate. Second, we identify a core as a connected 
region whose density is larger than = 10^ cm~^. Third, we subtract a center-of-mass 
velocity of the core from its velocity field. The three-dimensional center-of-velocity of the 
core ranges from 0.2 km sec~^ to 2.0 km sec~^ Fourth, we take the average of the density 
and radial velocity at each spherical shell with the width of a computational cell. At this 
point we have to mention the effect of the radial averaging that wc have taken in the above. 
The averaging certainly smoothed out density and velocity variations at an outer layer af- 
fected by, for example, a shock. So the line profiles observed towards a specific direction 
may show much stronger variation than the line profiles that are obtained after taking the 



-6- 



radial averaging. 

In Figure 1, we plot the density and velocity profiles as a function of the radial distance 
from the center of the core. The two digit numbers are the elapsed time in units of 40,000 
years after turning on self-gravity. We can clearly see from density profiles that the central 
density of the core is increasing as time goes on. The power index of density profiles in the 
outer region is increasing from shallower than -2 at an earlier time to steeper than -2 at a 
later time. At a given radial distance, the density and velocity dispersions are much smaller 
than the radially averaged density and velocity, respectively. 

The radial velocity profiles are quite interesting. In the central part of the core, most of 
the profiles have negative but very small velocity, which means that the central part of core 
is slowly collapsing. In the outer part of the core, the profiles have both positive and negative 
signs and their amplitudes are larger than the thermal sound speed. The large-amplitude 
velocity profiles are due to not only the large-scale turbulent flow nearby the core but also 
accretion shocks onto the core. Even though we do not show radial density and velocity 
profiles from a simulation with a weak initial field strength, (5 — 1.0, the large- amplitude 
velocity profiles in the outer part of a core are also found. However, infall velocities in the 
central part of the core have a bit larger (but still have a subsonic speed) than ones of the 
case with /3 = 0.1. 

Figure 2 shows the comparison between dust and gas temperatures calculated based 
on the dust continuum radiative transfer and the balance between heating and cooling of 
gas, respectively. The gas temperature is decoupled from the dust temperature at densities 
< 5 X 10^ cm~^, where gas- grain collisions is not dominant, and has peaks at the surface of 
the model core due to the photoelectric heating. 

The evolution of the HCO+ radial abundance calculated based on the density and tem- 
perature evolution is presented in Figure 3. As expected, at small radii, HCO'*' is signifi- 
cantly depleted from the gas phase as density grows because CO becomes frozen-out onto 
grain surfaces (Lee et al. 2003, 2004 and references therein) in such dense and cold con- 
ditions. The abundance drops at the core surface due to the dissociative recombination of 
HCO^ with electrons. As described in §2.2, the chemical evolution is calculated at each 
grid point without following each gas parcel (that is, assuming the quasi-static evolution) 
due to the following reasons. In the inner region of the core, the static evolution is a good 
approximation since the velocity is much smaller than the sound speed. In the outer part 
of the core, however, densities are low, and thus the chemical time scale is too long to be 
affected by the movement of gas parcels. For example, an infalling gas parcel with the infall 
velocity of 0.5 km s^^ at the surface moves inward only ~0.02 pc for 40,000 years, and the 
densities of two positions are not very different. This assumption is also well supported in 
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Figure 3. The abundance of HCO+ varies significantly at small radii, but at radii greater 
than 0.1 pc, where velocities are possibly greater than the sound speed, the abundance does 
not vary much with time. 

Figure 4 presents the evolution of the HCO+ 3—2, 4—3, and C^^O 3—2 line profiles 
toward the core center. The two digit number in each panel represents the elapse time in 
units of 40,000 years after turning on the self-gravity. The first upper left panel shows line 
profiles at 1 x 10^ (25 x 40, 000) years and the time of the second panel is increased by 
8 X 10'^ (2 X 40, 000) years from the first one, and so on. Although the density grows with 
time, the C^^^O 3—2 line becomes weaker since CO is depleted from gas due to its frcczc-out 
onto grain surfaces. The most important trend seen in this study is that the asymmetry 
of the HCO+ 3—2 line profile mostly depends on the velocity structure at the outer part of 
the core. In the comparison of hue profiles with the velocity structures (Figure lb), if the 
velocity at large radii shows infall motion(e.g., at the time of 37), the HCO"^ line has the 
infall asymmetry with a blue stronger peak. On the other hand, the line profile represents the 
reversed asymmetry with a red stronger peak when the velocity at large radii shows outward 
motion (e.g., at the times of 25, 31, and 49). This result is caused by the combination of 
the abundance and velocity profiles; HCO"'" is depleted at small radii so that the HCO"'" 3—2 
lines effectively form at rather large radii in the core, where velocities are dominated by the 
turbulent flow and accretion shocks. The kinetic temperature also increases outward. As a 
result, the excitation temperature of HCO+ 3—2 peaks around 0.1 pc. Even the broad fine 
feature seen at time 43, which may be interpreted as an infall at the core center in actual 
observations, is caused by the motion at large radii of the core rather than that occurring 
at small radii. 

In the comparison of the HCO''" 4—3 line (gray) with the 3—2 line (black sohd) in 
Figure 4, the 3—2 line shows the reversed asymmetry while the blue asymmetry appears in 
the HCO"*" 4—3 line at time 25 and 29 (the opposite feature at time 35 and 47) although the 
difference between blue and red peak strengths is not significant. This varying asymmetry 
in different transitions of a molecule has been observed (e.g., LI 157 in Cregersen et al. 1997; 
B68 in Maret et al. 2007). In addition, the relative strength between two peaks of the HCO"'" 
3—2 line is reduced in the 4—3 line at the times of 31 and 39. These trends are featured 
since the two transitions trace different radii of the core, for example, at time 25, the 4—3 
line traces a deeper region, which has inward motion, while the 3—2 line traces the less dense 
outer region, which shows outward motion. 

We assume that our model core of a radius 0.25 pc is at a distance 250 pc. Its angular 
size is about 200" . Since physical conditions of the core at around 0.1 pc (~ 80" ) mainly 
contribute to the line formation, a telescope beam size that is not larger than the angular 
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size of the core does not significantly affect tfie line intensity and shape shown in Figure 4. 

4. Summciry and Discussion 

We traced, through a high-resolution numerical experiment, the evolution of a core form- 
ing in a turbulent molecular cloud that is in a Jeans unstable and magnetically supercritical 
state. The core is slowly contracting at the central part with a subthermal speed. Due to 
the influence of accretion shocks and a nearby turbulent flow, velocity flelds in the outer 
region of the core can easily exceed the thermal speed, and show not only the contracting 
but also sometimes expanding motions. Based on the evolutionary profiles of density and 
gas and dust temperatures, we calculated the evolution of the HCO^ abundance. Once the 
central density of the core is larger than 10^ cm^'^, the HCO+ abundance profiles peak at 
outer regions because of the depletion of HCO+ at the central part of the core (see Figure 3). 
We also calculated the evolution of HCO"*" 3—2, 4—3, and C^^O 3—2 fine profiles that are 
coupled to the density, kinetic temperature, and abundance structures. The most interesting 
result in this work is that the HCO+ line profiles are dominantly affected by the velocity 
profiles at the outer regions. This is the consequence of both the HCO+ abundance and 
velocity profiles that have peaks at the outer regions. 

We showed that the molecular line profiles of a core, which results from the coupled 
velocity and chemical structures, vary with time. To understand the formation process of 
cores, survey observations of many cores in different evolutionary states are crucial. Mapping 
observations of individual cores with a high resolution enables us to know the complex 
velocity structures of the cores. The radial variation of the velocity field as seen in B68 
may be an intrinsic part of star formation, and the blue asymmetry is not always a direct 
signature of collapse in starless dense cores (Keto et al. 2006). 

We arc very grateful to Ted Bergin for valuable comments. This work was supported 
by the Korea Science and Engineering Foundation under a cooperative agreement with the 
Astrophysical Research Canter for the Structure and Evolution of the Cosmos. 
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Fig. 1. — The evolution of density and velocity profiles. Two digit numbers represent time 
in units of 40,000 years. 
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Fig. 2. — The evolution of the gas (thin lines) and dust (thick lines) temperature profiles. 
The two digit numbers represent time in units of 40,000 years. 
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Fig. 3. — The evolution of the HCO"'" abundance profiles. The two digit numbers represent 
time in units of 40,000 years. 
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Fig. 4. — The evolution of molecular line profiles. The black and gray sohd histograms 
represent the HCO"*" 3—2 and 4—3 line profiles, respectively. The dotted lines show the 
C^^^O 3-2 hnes. The beam sizes are 26" for HCO+ 3-2 and C^^O 3-2 and 20" for HCO+ 
4—3, respectively. The angular size of our model core with a radius of 0.25 pc at a distance 
of 250 pc is around 200" . The two digit numbers represent time in units of 40,000 years. 



